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Abstract 

An  analysis  has  been  developed  to  predict  transient  aeroelastic  response  of  gimballed 
tiltrotors  during  shipboard  engage/disengage  operations.  The  blades  are  modeled  as 
slender  elastic  beams  attached  to  a  gimballed  rotor  hub  undergoing  flap  and  lag  bending, 
elastic  twist,  and  axial  deflection.  The  gimbal  restraint  is  simulated  using  a  conditional 
rotational  spring.  The  rotor  equations  of  motion  are  formulated  using  Hamilton’s 
principle  and  spatially  discretized  using  the  finite  element  method.  The  discretized 
equations  of  motion  are  integrated  in  time  for  a  specified  rotor  speed  profile.  Blade 
element  theory  is  used  to  calculate  quasi-steady  loads  in  linear  and  nonlinear  regimes. 
Studies  for  a  l/5“*-size  aeroelastically  scaled  tiltrotor  are  conducted  to  validate  the  analysis 
and  investigate  transient  response  and  loads.  Good  correlation  existed  between  the 
experimental  data  and  the  prediction  of  blade  flap  and  lag  moments  in  a  hover  condition. 
Blade  bending  moment  and  hub  moment  predictions  indicated  that  gimbal  restraint 
impacts  could  induce  high  transient  loads  on  the  rotor  blades  and  hub.  A  new  method  for 
simulating  ship  airwake  flow  fields  is  presented.  Ship  airwake  flow  fields  are  inherently 
unsteady  and  very  difficult  to  predict  with  numerical  tools.  The  method  presented  in  this 
paper  is  called  the  Nonlinear  Disturbance  Equations  (NLDE)  and  is  fourth  order  accurate 
in  both  space  and  time.  In  this  method  the  steady  state  flow  field  is  solved  first,  and  then 
the  unsteady  fluctuations  are  solved.  Steady  and  unsteady  results  are  presented  for  a 
generic  frigate  shape.  The  parallel  method  is  a  necessity  for  ship  airwake  problems  and 
MPI  is  used  in  the  NLDE  solver.  Comparison  of  the  parallel  performance  on  various 
computers  is  also  presented. 
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TRANSIENT  RESPONSE  ANALYSIS  OF  GIMBALLED  TILTROTORS 
DURING  ENGAGE  AND  DISENGAGE  OPERATIONS 


Abstract 

An  analysis  has  been  developed  to  predict  transient  aeroelastic  response  of  gimballed 
tiltrotors  during  shipboard  engage/disengage  operations.  The  rotor  blades  are  modeled 
with  slender  elastic  beams  undergoing  flap  bending,  lag  bending,  elastic  twist,  and  axial 
deflection  attached  to  a  gimballed  rotor  hub.  The  gimbal  restraint  is  simulated  using  a 
conditional  rotational  spring.  The  rotor  equations  of  motion  are  formulated  using 
Hamilton’s  principle  and  spatially  discretized  using  the  finite  element  method.  The 
discretized  rotor  equations  of  motion  are  integrated  in  time  for  a  specifled  rotor  speed  run¬ 
up  or  run-down  profile.  Blade  element  theory  is  used  to  calculate  quasi-steady  loads  in 
linear  and  nonlinear  regimes.  Studies  for  a  l/5^'*-size  aeroelastically  scaled  tiltrotor  model 
are  conducted  to  validate  the  analysis  and  investigate  transient  response  and  loads  of  the 
gimballed  rotor.  Good  correlation  existed  between  the  experimental  data  and  the 
prediction  of  blade  flap  and  lag  moments  in  a  hover  condition.  Blade  bending  moment  and 
hub  moment  predictions  indicated  that  gimbal  restraint  impacts  could  induce  high 
transient  loads  on  the  rotor  blades  and  hub.  For  jx  =  0.096,  the  peak  flap  and  lag  bending 
moment  in  the  blade  flexure  are  130%  and  310%  above  loads  for  a  Ig  hover  condition. 

Introduction 

Unique  and  often  hazardous  conditions  are  encountered  when  rotorcraft  are  operated 
from  ship  based  platforms.  One  of  the  more  troublesome  rotorcraft/ship  interface  problems  can 
occur  during  the  engagement  and  disengagement  of  the  rotor  system  while  the  aircraft  is  on  the 
flight  deck.  Excessive  aeroelastic  flapping  due  to  high  wind  over  deck  conditions  can  occur  at 
low  rotor  speeds  when  the  centrifugal  stiffening  is  low.  For  articulated,  the  rotor  blade  can 
deflect  several  feet  and  contact  the  fuselage  of  the  helicopter.  For  gimballed  tiltrotors,  the  entire 
rotor  assembly  can  tilt  and  contact  the  gimbal  restraint,  possibly  resulting  in  high  blade  loads. 

Excessive  aeroelastic  flapping  of  rotor  blades  at  low  rotational  speeds,  or  “blade  sailing”, 
has  been  the  subject  of  both  analytical  and  experimental  investigations  over  the  past  several 
years*'**^.  In  the  United  Kingdom,  Newman  developed  an  analysis  for  the  transient  response  of 
rotor  blades  during  engagement  and  disengagement  operations.  A  set  of  unique  engage  and 
disengage  experimental  tests  were  conducted  by  using  scaled  models  in  a  wind  tunnel*' .  The 
analytical  model  was  applied  to  teetering  and  articulated  rotors.  At  Penn  State,  Geyer  and  Smith 
also  developed  a  transient  aeroelastic  response  analysis  for  helicopters  with  articulated  or 
hingeless  rotor  systems^'*.  The  finite  element  method  was  used  to  model  the  blade  as  an  elastic 
beam  undergoing  deflections  in  flap  bending  and  torsion.  Blade  element  theory  was  used  to 
calculate  the  quasi-steady  or  unsteady  aerodynamic  loads  in  linear  and  nonlinear  regimes. 

Present  tiltrotor  configurations,  such  as  the  V-22  Osprey,  feature  the  use  of  a  gimballed 
rotor  system.  A  gimballed  rotor  has  three  or  more  blades  attached  to  the  hub  without  flap  and 
lag  hinges.  The  hub  is  attached  to  the  rotor  shaft  by  a  universal  joint,  or  gimbal.  Other  unique 
rotor  characteristics  include  a  gimbal  spring  and  a  gimbal  restraint  as  shown  in  Figure  1.  In 
addition,  gimballed  rotor  blades  are  typically  short  and  highly  twisted. 


Hub  Undefiected  Hub  Deflected 


Figure  1  -  Gimballed  Tiltrotor  Hub  Schematic 

Compared  to  articulated  or  hingeless  rotor  systems,  the  gimballed  rotor  has  some  unique 
structural  and  aerodynamic  characteristics  related  to  shipboard  engagement  and  disengagement 
operations.  The  inflow  distribution  associated  with  the  ship  airwake  on  all  blades  contributes  to 
gimbal  flapping  moments.  The  motion  of  all  blades  is  coupled  through  the  rotor  hub.  In 
addition,  the  highly  twisted  blades  introduce  significant  structural  flap-lag  couplings  and  can 
result  in  high  angles  of  attack  at  low  rotor  speeds. 

An  analytical  model  for  steady-state  aeroelastic  response  and  stability  of  gimballed  tiltrotors 
was  developed  in  Ref.  11.  The  rotor  model  included  coupled  flap-lag  bending  modes  and  torsion 
degrees  of  freedom.  The  gimbal  tilt  was  described  by  using  the  gimbal  degrees  of  freedom,  Pgc 
and  Pgs,  which  are  the  pitch  and  roll  of  the  rotor  disk  in  nonrotating  frame.  In  1993,  Nixon 
developed  a  comprehensive  analysis  for  aeroelastic  response  and  stability  of  tiltrotor  with 
composite  rotor  blades^^.  Multi-blade  rotor  dynamics  were  studied  in  Ref.  13.  A  model  was 
developed  by  using  a  set  of  nonlinear  coupled  flap-lag-torsion  equations  of  motion  for  a  two- 
bladed  teetering  rotor.  The  sensitivity  of  aeroelastic  stability  boundaries  to  structural  and 
mechanical  coupling  was  explored.  For  these  analyses,  the  steady  periodic  response  of  the  rotor 
could  be  obtained  by  considering  the  response  of  an  equivalent  single  blade.  For  transient 
engagement  and  disengagement  analyses,  the  motion  of  blades  is  non-periodic.  In  addition, 
gimbal  restraint  impacts  make  the  individual  blades  respond  differently  and  the  inflow 
distribution  on  all  blades  contributes  to  gimbal  motion.  The  entire  multi-bladed  gimballed  rotor 
has  to  be  modeled  to  analyze  the  transient  response  during  engage  and  disengage  operations. 

Objectives 


The  overall  goal  of  this  research  is  to  systematically  investigate  the  transient  response  of 
gimballed  tiltrotors  during  shipboard  engage/disengage  operations  by 

1.  Developing  a  comprehensive  gimballed  rotor  model  based  on  an  elastic  finite  element  model. 

2.  Determining  transient  gimbal  motion,  blade  loads  and  hub  loads. 

3.  Conducting  parametric  studies  of  rotor  properties. 

A  separate  goal  of  this  research  project  is  to  develop  computational  techniques  to 
determine  the  airwake  environment  around  complex  ship  shapes.  The  results  of  this  research  are 
detailed  in  a  report  attached  to  this  one.  They  can  also  be  accessed  online  at 
http://www.personal.psu.edU/facultv/l/n/lnl/frame2.html#papers. 


Approach 


To  accurately  predict  the  blade  and  hub  loads  associated  with  a  gimbal  restraint  impact 
during  engage  and  disengage  operations;  an  elastic  blade  model  is  formulated.  The  following 
section  presents  the  methodology  and  results  from  this  analysis.  A  fixed  rotor  support  was 
considered  in  the  present  research  (i.e.  rigid  wing/pylon).  The  multi-bladed  gimballed  rotor  is 
modeled  as  several  slender  elastic  beams  attached  to  a  hub  and  undergoing  flap  bending,  lag 
bending,  elastic  twist,  and  axial  deflection.  The  blades  are  discretized  into  a  number  of  beam 
elements  using  the  finite  element  method.  Each  beam  element  consists  of  fifteen  degrees  of 
freedom  .  The  motion  of  gimbal  is  expressed  by  using  two  degrees  of  freedom  in  the  rotating 
system.  The  partial  differential  equations  governing  the  rotor  motion  are  derived  using 
Hamilton’s  Principle 


Sn  =  J {ST-SU-SW)dt  =  0 


(1) 


where  ST  is  the  variation  of  kinetic  energy,  SU  is  the  variation  of  strain  energy,  and  SW  is  the 
virtual  work  due  to  external  forces.  The  contributions  to  these  energy  expressions  may  be 
summed  as 


su=^su:  +SU^,  +SU^,  +SU^,  (2) 

m=l 

(3) 

m=\ 

Sw  =  '^dw,”'  (4) 

m=l 


Where  b  refers  to  the  blade,  gs  to  the  gimbal  spring,  gr  to  the  gimbal  restraint,  and  pi  to  the  pitch 
link.  Assuming  isotropic  elastic  properties  of  blades,  the  variation  of  blade  strain  energy  of  the 
blade  can  be  expressed  as 


=  f  0  )dTjdC  dr 


(5) 


where  is  the  axial  strain,  and  £xn  and  are  engineering  shear  strains,  E  is  the  axial  modulus 
of  elasticity,  and  G  is  the  shear  modulus.  The  variation  of  gimbal  spring  strain  energy  is  given 
by 


(6) 

The  degrees  of  freedom  in  the  rotating  system  {fi  and  are  used  to  express  the  motion  of 
gimbal.  The  transformation  matrix  is  given  by 


(7) 


|/^GC  1  [cos 

ly^GsJ 


-sin^  [p 
cosy/  1  ^ 


Where  y/is  the  azimuth  angle  of  rotating  frame,  P  and  ^are  the  gimbal  tilt  degrees  of  freedom  in 
the  rotating  system.  The  gimbal  degrees  of  freedom  in  nonrotating  and  rotating  coordinate 
systems  are  shown  in  Figure  2. 


Figure  2  -  Gimbal  Degrees  of  Freedom 


The  variation  of  gimbal  spring  strain  energy  in  the  rotating  system  is  much  like  Eqn.  6 

SU^^=KppSp+Kp(l>8(t>  (8) 


This  analysis  models  the  gimbal  restraint  interaction  with  the  rotor  using  a  conditional  rotational 
spring  in  nonrotating  frame.  The  rotational  spring  stiffness  is  zero  for  tilt  angles  less  than  the 
gimbal  restraint  angle  and  large  enough  to  restrict  further  tilt  to  less  than  0.1°  for  tilt  angles 
larger  than  the  gimbal  restraint  angle.  The  variation  of  strain  energy  due  to  the  gimbal  restraint 
is  expressed  as 


=  (9) 

The  maximum  tilt  angle,  P^ax,  is  given  by 

Anax  =  {P  y/-y>sm  yf)cos  y/^^  +  {Ps\ny/  +  (l)cosy/)siny/^^  (10) 

Eqn.  9  can  be  expressed  in  the  rotating  system  as: 

5U^  =  [A^pdp+B'(l)d(l>+AB(l)5p+ABP6(l)-AP^6p-Bp^S<l))  (1 1) 


where 


A  =  cos  cos  y/^  +  sin  yr  sin  y/^^ 

B  =  -s\x\yf  cos  y/^^  +  cos  yr  sin  y/^^ 


(12) 


where  y/max  is  the  azimuth  angle  where  the  maximum  gimbal  tilt  angle  is  located  and  is  given  by 


(13) 


fisini/r  +  tpcosyr^ 

Pcosyz-^smyr  ^ 

The  virtual  work  can  be  expressed  as  a  summation  of  contributions  from  gravitational  forces  and 
aerodynamic  forces 


^<^max=tan 


SW  =  SWg+SW^p  (14) 

Gravitational  forces  are  required  in  this  analysis  due  to  their  effect  on  blade  deflection  over  the 
low  rotor  speed  regions.  The  work  performed  on  the  blade  due  to  gravitational  forces  is 
expressed  as 


SW(;=-{  mgdwdr  (15) 

J  0 

where  g  is  the  acceleration  due  to  gravity.  The  virtual  work  performed  by  the  aerodynamic 
forces  is  given  by 


=^’^^{l^5v  +  l!lSw  +  M^5$^dr  (16) 

The  unsteady  aerodynamic  model  developed  by  Leishman  and  Beddoes  is  used  for  blade  element 
airloads.  Significant  documentation  on  the  development  of  this  theory  can  be  found  in  Refs.  15 
and  16.  The  model  is  divided  into  three  flow  regions,  attached  flow  unsteady  aerodynamics, 
nonlinear  aerodynamics  (trailing  edge  separation),  and  dynamic  stall.  The  unsteady  attached 
flow  model  is  based  on  an  indicial  response  formulation  and  predicts  both  noncirculatory  and 
circulatory  airloads.  A  detailed  description  of  the  incorporation  of  Leishman  and  Beddoes’ 
unsteady  aerodynamic  model  into  the  present  analysis  is  given  in  Ref.  6.  During  engage  and 
disengage  operations,  the  rotor  speed,  Q,  is  very  low.  The  blades  can  experience  large  reverse 
flow  regions  and  very  high  angles  of  attack  compared  to  operation  at  normal  rotor  speed.  The 
aerodynamic  model  used  in  this  analysis  must  be  modeled  carefully  for  this  situation.  In  order  to 
simulate  the  aerodynamic  loads  at  very  high  angles  of  attack  (|«|  >  20°),  a  semi-empirical  quasi¬ 
steady  model  was  used  to  improve  the  quasi-steady  aerodynamic  force  prediction*^.  For  the 
angle  of  attack  of  section  higher  then  20°,  the  section  lift,  drag,  and  pitching  moment  coefficients 
are  expressed  as 


c,  =1.15(l  +  M)‘’®sin2(a-a^) 

Crf  =1.05-h0.5M-1.04cos2a  (17) 

=-0.55sin(ar-10) 

where  oris  the  section  angle  of  attack,  c^o  is  the  zero  lift  angle,  and  M  is  the  Mach  number.  A 
quasi-steady  aerodynamic  model  is  used  in  this  analysis.  A  linear  gust  model  is  used  to  simulate 
the  airflow  over  a  ship  flight  deck^.  Detailed  descriptions  of  the  quasi-steady  aerodynamic 
modeling  and  the  shipboard  aerodynamic  environment  are  given  in  Ref  18  and  Ref  6 
respectively. 


The  variation  of  kinetic  energy  is  given  by 


ST=j’Jlp.V^Vd^d(dr 


(18) 


where  Ps  is  the  mass  density  of  the  blade  and  V  is  the  velocity  vector  at  any  point  on  the  blade 
section.  For  engage  and  disengage  operations,  the  rotating  speed  of  the  rotor  is  a  function  of 
time.  The  angular  acceleration  of  the  rotor,  Q ,  was  also  included  in  the  kinetic  energy. 

An  illustration  of  the  blade  finite  element  discretization  is  shown  in  Figure  3.  Each  of  the 
flexible  elements  consists  of  fifteen  degrees  of  freedom.  These  degrees  of  freedom  are 
distributed  over  five  element  nodes,  two  external  and  three  internal,  which  describe  the  elemental 

flap  deflection  (w,  w’),  lag  deflection  (v,  v),  twist  (^),  and  axial  deflection  (u).  Deformations 
within  the  element  are  described  by  shape  functions.  A  more  detailed  description  of  these 
degrees  of  freedom  and  shape  functions  is  given  in  Ref.  18.  Note  that  the  pitch  bearing  is  located 
at  the  outboard  part  of  the  flexure. 


Pitch  Bearing 


Figure  3  -  Finite  Element  Model  of  Gimballed  Rotor 

Figure  4  shows  the  global  degrees  of  freedom  for  the  gimballed  rotor.  The  first  two 
degrees  of  freedom  describes  the  gimbal  tilt.  The  last  degrees  of  freedom  at  the  left-hand  node 
of  the  second  element  of  each  blade  describe  the  rotation  of  the  pitch  bearing.  The  second 
element  is  the  flexure  of  the  blade.  For  ten  elements  in  each  blade,  there  are  275  DOF  with  30 
elements  for  the  entire  rotor. 

The  energy  expressions  are  spatially  discretized  by  substituting  the  shape  functions  into 
the  elemental  virtual  energy  expressions.  The  resulting  expressions  are  integrated  in  space  using 
a  six  point  Gaussian  Quadrature  method.  The  virtual  energy  expression  for  multi-bladed  tiltrotor 
in  terms  of  the  elemental  matrices  and  load  vector  becomes 


Figure  4  -  Global  Degrees  of  Freedom  and  Finite  Element  Model  Descriptions 


=  J i  E  -jqr '  (Mrqr + crqr + Kfqr  -  f”  )dt + 
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where  Kg  is  the  stiffness  of  matrix  due  to  gimbal  restraint  and  is  expressed  as 


Kp  0 

0  K„ 


(20) 


Nb  is  the  number  of  blades  for  tiltrotor  and  N  is  the  total  number  of  spatial  finite  elements  in  one 
blade.  The  element  nodal  displacement  vector,  qf ,  is  defined  as 

q”'={<,  «J  M,  V,  v'  Vj  v'  W,  VI''  Wj  w'  02  03]^  (21) 

Where  the  subscript  i  refers  to  the  element  of  the  m'*  blade.  For  a  gimballed  tiltrotor,  the  blade 
is  fixed  to  the  hub  at  the  blade  root.  The  hub  is  attached  to  the  rotor  shaft  by  a  gimbal.  The 
motion  of  the  hub  is  expressed  using  two  degrees  of  freedom  in  this  analysis,  P  and  which 
rotate  with  the  rotor.  The  motion  at  the  root  of  the  blades  can  be  determined  by  six  degrees  of 

freedom, and  For  the  gimballed  rotor  boundary  conditions,  the  axial, 
vertical,  and  chordwise  displacements  response  quantity,  u,  v,  and  w,  at  each  blade  root  are 
specified  to  be  zero.  The  lag  angle,  v’,  is  specified  to  zero  as  well.  Only  the  root  slope,  w,""  and 

twist,  are  not  restrained  and  can  be  expressed  as  a  function  of  two  gimbal  degrees  of 
freedom  in  the  rotating  frame,  y^and 


w["'  =ficos 


(22) 


(23) 


where  the  subscript  1  indicates  the  left-hand  node  of  the  first  blade  finite  element.  In  this 
analysis,  the  node  number  order  is  defined  from  the  root  to  the  tip  of  the  blade.  The  azimuth 
angle  of  each  blade,  y/m,  in  the  rotating  frame  can  be  expressed  as 

(24) 

Eqns.  22  and  23  are  then  substituted  into  Eqn.  19.  The  stiffness  matrix,  mass  matrix,  damping 
matrix  and  force  vector  for  each  element  are  then  assembled  to  form  their  global  counterparts. 
The  variation  of  the  total  energy  of  the  multi-bladed  gimballed  rotor  system  becomes 

<5n  =  f'"  (Mq  +  Cq  +  Kq  -  F)  rit  (25) 


Finally,  the  discretized  equation  of  rotor  motion  is  obtained  as 

Mq+Cq+Kq=F 


(26) 


where  M,  C,  K,  and  F  are  mass  matrix,  damping  matrix,  stiffness  matrix,  and  force  vector  of  the 
entire  rotor.  The  mass  matrix,  damping  matrix  and  stiffness  matrix  include  the  linear  structural, 
inertial  and  aerodynamic  contributions.  The  force  vector  includes  the  constant  structural,  inertial 
and  aerodynamic  contributions  as  well  as  the  nonlinear  structural  and  aerodynamic  contributions. 

During  gimbal  restraint  impacts,  it  is  important  to  calculate  blade  and  hub  loads  for 
design  limits.  The  flap  bending  moment,  lag  bending  moment,  and  twisting  moment  in  the 
deformed  coordinate  system  is  defined  by’® 
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where  A,  Aca,  AK/,  E/,  B2,  Ci,  and  C2  are  section  integrals.  The  pitch  angle,  0,  is  defined  as 


0  =  00+  0„  cos  Wm  +  cos  ¥m+0yr-^)-  k^w,”' 


(28) 


Elastic  Blade  Model  Results 


Validation 

Experimental  data  for  blade  bending  moments  and  twisting  moments  in  a  hover  condition 
of  a  l/5**''scale  model  of  the  V-22  gimballed  rotor  was  used  to  validate  the  analysis.  The  test  was 
conducted  for  a  hover  condition  with  varying  collective  pitch  and  zero  cyclic  pitch.  The  rotor 
speed  is  875  RPM.  The  flap  bending  moment,  lag  bending  moment,  and  twisting  moment  at  six 
blade  sections  was  measured.  A  comparison  of  these  experimental  results  with  the  predictions  of 
the  present  analysis  is  shown  in  Figure  8a  and  8b.  The  flap  and  lag  moments  predicted  by  this 
analysis  show  good  correlation  with  the  experimental  test  data  at  moderate  collective  pitch.  It 
should  also  be  noted  that  no  nonlinear  aerodynamics  was  used  for  these  calculations.  Note  that 
the  element  between  section  O.llR  and  0.157?  represents  the  coning-flexure  of  the  blade. 
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Figure  5a  -  Flap  Moment  Comparison 
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Figure  6a  -Lag  Moment  Comparison 


Transient  Response  Analysis 

The  blade  transient  response  is  calculated  by  integrating  the  discretized  rotor  equations  of 
motion  for  a  specified  rotor  speed  run-up  profile.  The  rotor  speed  and  azimuth  run-up  profiles 
are  prescribed  using  experimental  data.  The  experimentally  determined  rotor  speed  run-up 
profile  (for  the  first  20  seconds)  are  expressed  as 


Q(r) =3.08x10^  xt*  -1.15x10^  xt^  +1.79x10-^  xt" 
-3.03xl0r^  xr^  -4.92xl0r^  xt"  -1-0.164 


In  an  effort  to  clearly  understand  engage/disengage  behavior,  an  investigation  was  focused  on  the 
motions  of  the  tiltrotor  gimbal  and  bending  moments  of  the  blades.  The  tiltrotor  operation  was 
simulated  for  a  starboard  wind-over-deck  of  ^  =  0.075.  A  25%  upflow  through  the  windward 
half  of  the  rotor  disk  and  a  25%  downflow  through  the  leeward  half  were  simulated.  The 
corresponding  full-scale  wind  speed  is  equivalent  to  the  35  knots.  The  rotor  collective  pitch  is 
set  at  0°  with  zero  cyclic  control  input.  Figure  9  shows  the  time-history  of  transient  response  of 
gimbal  tilt,  P  and  and  the  maximum  tilt,  P^ax  for  the  run-up  operation.  The  azimuth  of  rotor  in 
graph  is  defined  as  the  azimuth  of  rotor  rotating  and  measured  from  downstream  (forward  flight) 
to  span  axis  of  the  blade  1.  The  maximum  tilt  indicates  the  impacts  between  rotor  and  gimbal 
restraint  during  the  run-up  operation.  The  flap  and  lag  bending  moment  time  histories  at  section 
r/7?=0.15  (inboard  section  of  flexure)  in  the  blade  1  are  shown  in  Figure  8.  It  can  be  seen  that 
impacts  between  rotor  and  gimbal  restraint  cause  substantial  increase  in  blade  bending  moments. 
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Figure  7  -Transient  Response 


Figure  8  -Bending  Moment  Response  in 
Flexure 


Wind  speed  effects 

Blade  moments  were  predicted  for  //  =  0.075  and  fx  =  0.096  starboard  winds  with  a  25% 
upflow  through  the  windward  half  of  the  rotor  disk  and  a  25%  downflow  through  the  leeward 
half  of  the  rotor  for  WRATS  model  rotor,  which  is  corresponding  to  35  and  50  knots  for  full- 
scale  V-22  rotor.  Rotor  collective  pitch  is  set  at  0°  with  zero  cyclic  control  input. 
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Figure  9  -  Comparison  of  Bending  Moment  at  section  0.15/?  For  Different  Wind  Speeds 

Figure  9  shows  the  comparison  of  blade  flap  and  lag  bending  moment  peak  predictions  at 
inboard  section  of  the  blade  flexure  (0.15/?)  during  run-up  operation.  It  shows  that  the  bending 
moment  peaks  during  run-up  operation  are  much  higher  than  in  hover  condition.  The  moment 
peaks  appear  when  gimbal  restraint  impacts  occur.  For  =  0.075,  the  flap  bending  moment  is 
30%  above  loads  for  IG  hover  condition;  the  lag  bending  moment  is  130%  above  loads  for  hover 
condition.  The  flap  bending  moment  peak  for  jx  =  0.096  is  1.3  times  bending  moment  peak  for  n 
=  0.075.  The  lag  bending  moment  peak  for  the  higher  wind  speed  is  1.8  times  moment  peak  for 
lower  wind  speed. 


Flap  Lag 


Summary  and  Conclusions 

An  analysis  has  been  developed  to  predict  transient  aeroelastic  gimbal  tiltrotor  response 
during  rotor  engage/disengage  sequences.  The  multi-blade  gimballed  rotor  is  modeled  as  several 
slender  elastic  beams  attached  to  a  hub,  which  undergo  flap  bending,  lag  bending,  elastic  twist, 
and  axial  deflections.  The  gimbal  restraint  is  simulated  using  a  conditional  rotational  spring. 


Gimbal  motion  is  expressed  by  using  two  degrees  of  freedom  in  the  rotating  frame.  The  rotor 
equations  of  motion  were  formulated  using  Hamilton’s  principle  and  spatially  discretized  using 
the  finite  element  method.  The  discretized  rotor  equations  of  motion  were  integrated  in  modal 
space  for  a  specified  rotor  speed  run-up  or  run-down  profile.  Blade  element  theory  was  used  to 
calculate  quasi-steady  loads  in  linear  and  nonlinear  regimes.  Parametric  studies  for  a  l/5‘''  size 
aeroelastically  scaled  V-22  tiltrotor  model  were  conducted  to  investigate  the  transient  response 
characteristics  and  effects  of  wind  speed  and  blade  flexure  stiffness  on  transient  rotor  motion  and 
blade  and  hub  loads.  Predictions  showed  that  for  =  0.075,  the  flap  bending  moment  is  30% 
above  loads  for  Ig  hover  condition;  the  lag  bending  moment  is  130%  above  loads  for  hover 
condition. 
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Abstract 

This  paper  presents  a  new  method  for  simulating  ship 
airwake  flow  fields.  These  flows  are  inherently  unsteady, 
and  very  difficult  to  predict.  The  method  presented 
(NLDE)  is  fourth-order  accurate  in  space  and  time.  In 
this  method  we  first  solve  for  the  steady  state  flow  field, 
then  we  solve  for  the  unsteady  fluctuations.  Steady  and 
unsteady  results  are  presented  for  a  generic  frigate  shape. 
The  parallel  method  is  a  necessity  for  ship  air  wake  prob¬ 
lems  and  MPI  is  used  in  the  NLDE  solver.  The  parallel 
performance  on  various  computers  is  compared  also. 

Introduction 

Sharp-edged  box-like  ship  super-structures  create  nu¬ 
merous  aerodynamic  and  fluid  dynamic  problems.  Un¬ 
steady  separated  flow  from  sharp  edges  (and  excessive 
ship  motions)  make  landing  helicopters  on  ships  a  very 
hazardous  operation.  In  addition,  the  strong  unsteady 
flows  can  cause  severe  rotor  blade  deformations.  There 
have  been  numerous  incidences  where  the  helicopter 
blades  have  actually  impacted  the  helicopter  fuselage, 
which  is  called  a  “tunnel  strike”  In  order  to  avoid 
this  and  other  engage/disengage  problems,  determining 
safe  operating  envelopes  is  very  costly  and  time  consum¬ 
ing.  On  the  other  hand,  many  numerical  simulation  at¬ 
tempts  of  this  flow  field  have  not  been  successful  due 
to  the  inherently  unsteady  nature  of  flow  and  the  low- 
speed  character  of  the  flow  (which  may  cause  numerical 
stiffness). 

Research  on  ship  airwakes  has  been  conducted  using 
several  different  approaches.^  One  of  the  sources  of  rele¬ 
vant  research  is  building  aerodynamics  which  shows  the 
general  features  of  flow  about  blunt  bodies  of  different 
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aspect  ratios.  The  simplest  model  of  a  ship,  admittedly 
rather  crude,  is  a  sharp  edged  blunt  body.  The  super¬ 
structure  of  most  modern  ships  is  very  complicated,  in¬ 
cluding  towers,  antennae,  radar  dishes,  exhaust  stacks, 
etc.  The  flow  around  these  obstacles  is  very  difficult  to 
predict. 

Geometrically  precise  studies  are  needed  and  have 
been  done  in  wind  tunnels.^’ There  have  also  been 
full  scale  tests  performed  by  the  US  Navy,^^  which  gives 
some  important  information  on  real  ship  airwakes.  Of 
course  it  is  difficult  to  perform  very  controlled  experi¬ 
ments  on  real  ships.  It  is  also  difficult  to  measure  the 
flow  field  accurately  in  the  harsh  ocean  environment  and 
in  the  presence  of  the  strong  electromagnetic  fields  on 
most  ships. 

Most  wind  tunnel  tests  include  measurements  made 
in  the  wake  of  a  model  ship  exposed  to  a  uniform  veloc¬ 
ity  profile  and  almost  zero  turbulence  level.  One  more 
realistic  test  was  conducted  at  NASA  Ames  in  the  “Ship¬ 
board  Simulator”  with  a  neutrally  buoyant  atmospheric 
boundary  layer. ^ 

Another  reference  for  simulations  is  that  by  NRC- 
CNRC.^^  A  wind  tunnel  investigation  of  the  characteris¬ 
tics  of  the  airwake  behind  a  model  of  a  generic  frigate  was 
conducted.  The  wind  tunnel  simulation  incorporated  a 
correctly-scaled  atmospheric  boundary  layer.  Measure¬ 
ments  of  streamwise  and  vertical  components  of  airwake 
velocity  were  made.  Time  average,  standard  deviations,* 
spectral  densities  and  time  correlations  are  presented  for 
both  velocity  components  for  various  position  in  the  air¬ 
wake. 

The  wind  tunnel  tests  have  to  suitably  scale  the  en¬ 
vironment  and  structure  to  model  size,  make  the  appro¬ 
priate  measurements  in  the  wake  of  the  model  and  then 
rescale  the  results  back  to  full  size.  All  these  experi¬ 
mental  tests  are  crucial  for  validating  numerical  models. 
Wind  tunnel  tests  can  be  quite  costly,  but  flow  measure¬ 
ments  on  real  Naval  ships  are  very  difficult  and  costly  to 
obtain. 

Figs,  1  and  2  show  a  frigate  and  an  LHA,  respec- 


tively.  These  are  very  different  ships,  and  their  air  wakes 
are  very  different  also.  The  frigates  typically  carry  one 
or  two  SH-2G  Seasprites  or  SH-60B  Seahawks.  On  the 
frigate  we  are  mainly  interested  in  studying  the  hangar 
deck  area  (aft  portion  of  the  ship),  and  the  separated 
flow  that  effects  this  region.  On  the  LHA,  helicopters 
can  land  on  many  different  locations  on  the  deck,  and 
each  of  these  can  experience  quite  different  flow  fields. 
The  LHA’s  can  carry  9  CH-53D  Sea  Stallions  or  12  CH- 
46D  Sea  Knight  helicopters,  and  6  AV-8B  Harriers.  The 
forward  portion  of  the  deck  is  primarily  influenced  by 
the  separated  flow  off  the  deck  edge.  Very  strong  vortex 
sheets  emanate  from  these  edges.  One  of  the  authors 
(Long)  spent  three  days  on  an  LHA  (U.S.S.  Saipan)  and 
helped  Kurt  Long  measure  ship  airwakes.  We  found  in 
some  cases  the  flow  velocity  ranged  from  40  knots  12  feet 
off  the  deck  to  zero  velocity  3  feet  off  the  deck.  In  the 
mid-section  of  the  ship  the  very  large  island  has  a  strong 
effect  on  the  flow  and  tunnels  the  flow  tangential  to  the 
island. 


Figure  1:  Oliver  Hazard  Perry  Class  Guided  Missile 
Frigate  (length=445  feet,  beam=45  feet) 

The  need  for  numerical  simulations  comes  from  the 
very  high  cost  of  determining  the  safe  operating  en¬ 
velopes  for  helicopters  in  a  ship  environment  (and  the 
huge  testing  backlog).  It  would  be  very  useful  to  have 
numerical  methods  that  could  accurately  simulate  ship 
airwakes.  There  have  been  other  attempts  at  numer¬ 
ically  simulating  ship  airwakes.  The  air  wake  about  a 
DD-963  ship  configuration  was  simulated  using  a  steady- 
state  flow  solver  based  on  the  3D  multi-zone,  thin-layer 
Navier-Stokes  method.^®  A  US  navy  destroyer,  DDG51 
was  chosen  to  validate  an  unsteady  inviscid  solver  with 
an  unstructured  grid  and  low-order  method. No 
method  to-date  has  been  entirely  satisfactory  for  pre¬ 
dicting  these  flow  fields. 

Flow  Nature  of  Ship  Airwake 


Figure  2:  Tarawa  Class  LHA  (length  =  820  feet,  beam 
=  132  feet) 

Simulation 

From  previous  studies,  it  has  been  shown  that  the 
key  features  of  ship  airwakes  are  (1)  a  low  Mach  num¬ 
ber  (about  0.05),  (2)  inherently  unsteady  flow,  and  (3) 
large  regions  of  separated  flow.  The  large  separated  re¬ 
gions  from  superstructure  sharp  edges  are  quite  difficult 
to  capture  accurately.  In  addition,  the  wind  conditions 
over  rough  seas  have  to  be  considered,  such  as,  the  at¬ 
mospheric  turbulent  boundary  layer  and  the  effect  of 
the  wind/ship  speed  ratio  on  the  turbulence  intensity. 
When  this  ratio  is  increased,  the  turbulence  intensity 
will  decrease  and  its  spectrum  will  shift  to  a  high  value 
in  the  streamwise  direction.  The  wind  direction  can  vary 
a  great  deal,  since  the  air  flow  can  impact  the  ship  at  any 
yaw  angle  (even  180  degrees).  The  complex  ship  geome¬ 
try  makes  unstructured  grid  solvers  and  parallel  comput¬ 
ers  very  attractive.  In  this  paper,  preliminary  attempts 
at  high  order  accurate  ship  airwake  predictions  have  been 
made  by  solving  a  steady  flow  field  with  a  well-developed 
CFD  method  (CFL3D^*)  and  a  perturbation  field  with  a 
high-order  method.  The  result  is  high-order- accurate  3D 
simulations.  We  will  only  show  structured  grid  results 
herein,  but  we  are  also  pursuing  higher-order  accurate 
unstructured  solvers  (e.g.  based  on  PUMA^). 

Nonlinear  Disturbance  Equations 
(NLDE) 

The  methodology  used  here  is  based  on  the  nonlinear 
disturbance  equations,  which  is  a  newly  developed  nu¬ 
merical  method, The  general  Navier-Stokes  equations 
in  a  Cartesian  coordinate  system  are: 

dt  ^  dx  ^  dy  ^  dz  dx  dy  dz 


(1) 


where  F,  G,  and  H  are  the  inviscid  terms  and  R,  S, 
E  are  the  viscous  terms.  The  results  presented  here  will 
all  be  inviscid.  The  flow  field  is  then  split  into  a  mean 
and  a  fluctuating  part: 
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Substitution  of  equation  (2)  into  (1)  and  rearranging 
results  in  the  nonlinear  disturbance  equations  (NLDE): 


^  ^  ^ 

dt  ^  dx  ^  dy  ^  dz 
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On  the  left  hand  side  of  the  NLDE  are  terms  related 
to  the  perturbation  properties  and  the  cross  terms  (lin¬ 
ear  and  nonlinear),  whereas  the  right  hand  side  contcdns 
strictly  mean  flow  terms. 

The  convective  fluxes  involving  the  perturbation  quan¬ 
tities  F\G*  and  are  given  as 
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The  mean  flow  source  term  Q  is  time  independent: 
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If  the  NLDE  is  time  averaged,  it  becomes  the 
Reynolds-averaged  Navier-Stokes  equation,  where  the 
Reynold’s  stresses  are  on  the  left  hand  side.  Thus,  for  a 
laminar  flow  <5  =  0. 

We  seek  a  solution  of  the  perturbation  variables  q* 
with  a  known  mean  flow  field  which  can  be  obtained 
from  existing  well-developed  CFD  codes  (e.g.  CFL3D,^® 
INS3D,  OVERFLOW,  PUMA,^  ...)  for  steady  flow. 
This  methodology  allows  us  to  use  the  most  effective 
algorithms  for  the  steady  and  unsteady  portions  of  field, 
respectively.  It  also  minimizes  round-off  error  since  we 
are  only  computing  perturbations.  We  can  even  use  dif¬ 
ferent  grids  for  the  steady  and  unsteady  solution.  More 
discussion  on  this  new  method  is  in  the  reference.^^ 


Characteristic  Boundary 
Conditions  for  NLDE 


The  boundary  conditions  for  the  NLDE  are  developed 
by  applying  Thompson’s  characteristic  method^^  to  the 
nonlinear  disturbance  equations.  Instead  of  using  the 
conservative  form  of  the  equations,  the  boundary  condi¬ 
tions  are  derived  based  on  the  nonconservative  form: 
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Here  VISp  and  VISvei  are  viscous  terms.  Substitu¬ 
tion  of  equation  (2)  into  (11)  -  (13)  and  rearranging  the 
results  gives  the  boundary  conditions  in  nonlinear  per¬ 
turbation  form: 
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Here  we  apply  a  characteristic  analysis  to  the  pertur¬ 
bation  variables  and  move  the  mean  flow  terms  to  the 
right  hand  side,  which  is  different  than  the  approach  em¬ 
ployed  by  others.  Now  considering  a  boundary  located  at 
X  =  xq  and  using  the  characteristic  analysis^^  to  modify 
the  hyperbolic  terms  of  Eqs.  (14)  -  (15)  corresponding 
to  waves  propagating  in  the  x  direction,  we  caxi  recast 
this  system  as: 


%  +  ^[^2  +  \(Lx  +  Ls)]  +  (po  +  P')(^  +  ^) 

Hvo  +  ^')^  +  (W'o  +^')%  =  +  ^  + 


r  j. 

W-57  + 


/  t\dp'  I  t.dp'  ,  nduo 

K+«)^  +  (»o  +  »)^=-(7p(-^  + 

5uo  ,  dwo  ,dpo  ,dpo  ^  ,dpo., 

+  +  ar» 

duo  ,  dvo  ,  dwo.,  dpo 


du'  1  r  \ 

+  2(po +7)^ -  Li)  +  (.0  + 1;  )-^  + 

,  au'  ,3wo  ,5uo  ,  ,5uo, 

9uo  ,  9«o  ,  ,  1  dpo 

ar>  +  ^^T7  a? 


+VJ5„  (19) 


dvo  dwo.  ,  dpo  dpo  ^P0\-if,’7\ 

■§5-  +  &->  +  <“»■&  +  ””%  +  ”»  9F>1 


-i  +  j(L.+ls)+7(p.+p')(-^  +  -57)  + 


+  I3  +  (^0 )-^  +  (tt'o +«»  )-^ + 
1  dp'  ,  ,5i;o  ,  ,dvo  ,dvo. 

,,  avo  ,  dvo  ,  1  9po 

+F7Sv 


dw'  ,  , 

-^  +  L4  +  (^^o  +  t^)^  +  («;o+«;)-^  + 

1  dp'  ,  ,dwo  ,  ,5wo  ,  ,dwQ 

..  3u;o  ,  dwo  ,  ,  1  9po 

“[(uo“^ - (-  vq-^ - h  zx;o-^)  -t - 

ox  oy  oz  Po  +  p'  dz 

^VIS^  (21) 

where  the  L^’s  are  the  amplitudes  of  characteristic 
waves  associated  with  each  characteristic  velocity  A^. 
These  velocities  are  given  by: 

Ai  =  (uo  '^u')-c  (22) 


A2  =  A3  =  A4  —  uq  +  u'  (23) 
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where  c  is  the  speed  of  sound: 
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Ai  and  A5  are  the  speed  of  acoustic  waves  moving  in 
the  negative  and  positive  x  direction;  A2  is  the  convection 


velocity  (the  speed  at  which  entropy  waves  will  travel) 
while  A3  and  A4  are  the  velocities  at  which  v  and  w  are 
advected  in  the  x  direction.  The  L,'s  are  given  by; 
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We  have  shown  here  the  characteristic  boundary  con¬ 
ditions  in  the  x  direction.  The  equations  for  the  y  and  z 
direction  are  derived  in  a  similar  manner.  For  edges  and 
corners  in  three  dimensional  situations  a  simple  exten¬ 
sion  is  used  by  combining  the  equations  for  the  various 
directions  into  one  approximate  equation. 

This  type  of  boundary  condition  treatment  also  al¬ 
lows  one  to  easily  introduce  a  disturbance  at  the  in¬ 
coming  boundary  by  deriving  an  expression  for  one  of 
the  incoming  characteristics  with  a  source  term.  Atmo¬ 
spheric  boundary  layer  conditions  can  be  incorporated 
in  this  manner  at  incoming  boundaries.  At  the  outflow 
boundaries,  the  boundary  conditions  are  essentially  non¬ 
reflecting.  The  ship  superstructure  and  ocean  surface  are 
both  treated  ais  hard  wall  boundary  conditions. 

Numerical  Method  and  Parallel 
Methodology 

The  NLDE  are  cast  in  a  generalized  coordinate  system 
and  solved  numerically  using  a  finite  difference  based 
scheme.  The  discretized  equations  are  solved  in  a  time 
accurate  manner  by  taking  advantage  of  computational 
aeroacoustics  (CAA)  methods.  The  spatial  flux  deriva¬ 
tives  are  calculated  using  seven  point  stencils  of  the 
fourth  order  optimized  Dispersion  Relation  Preserving 
(DRP)  scheme  of  Tam  and  Webb.^®  The  time  integra¬ 
tion  is  a  fourth  order  accurate  Runge-Kutta  method. 

Efficient  computing  performance  is  achieved  by  using 
a  three  dimensional  domain  decomposition  strategy.  The 
code  is  written  in  Fortran  77  plus  Message  Passing  In¬ 
terface  (MPI)^^  and  is  scalable  in  three  dimensions.  As 
mentioned  early,  the  ship  geometry  is  very  complicated, 
even  for  a  generic  frigate  test  model.  This  makes  multi¬ 
block  grid  simulations  and  domain  decomposition  very 
difficult.  In  order  to  make  the  code  scalable  and  flexible, 
a  three  dimensional  single-block  grid  is  used.  The  whole 


computational  domain  is  divided  into  many  three  di¬ 
mensional  zones.  The  grid  points  are  evenly  distributed 
across  each  processor. 

The  NLDE  solver  is  implemented  port  ably  on  paral¬ 
lel  computers,  such  as,  the  IBM  SP2  (e.g.  Penn  State, 
Npaci,  MHPCC),  SGI  Power  Challenge  and  Pentium 
II  Cluster.  A  comparison  of  code  performance  for  the 
ship  airwake  run  on  various  machines  is  shown  in  Fig.  3. 
While  a  24-processor  IBM  SP2  is  8.4  times  faster  than  8 
Pentium  IPs  networked  together,  the  SP2  costs  roughly 
14  times  more  than  the  PC  cluster.  Fig.  4  gives  the 
wall  clock  time  for  a  ship  air  wake  case  with  1.86  mil¬ 
lion  grid  points  using  various  number  of  processors.  A 
64- processor  SP2  is  roughly  2.6  time  faster  than  a  16- 
processor  SP2  (when  problem  size  is  kept  fixed). 


PC  Cluster 


1.  On  8  Pentium  tt  (266  MHz)  cluster  - 

with  1 00  Mb/sec  network  : 

256  MB  each  node  “ 

2.  On  0  R10000  nodes  (1 95  MHz) 

4  GB  shared  memory 

3.  On  24  Power2  thin  nodes  (66MHz)  - 

128  MB  each  node  : 

4.  On  24  Powef2  Super  Chips 
(160  MHz,  thin  node) 

512  MB  each  node  : 

5.  On  24  Power2  Super  Chips  j; 

(1 60  MHz,  thin  node) 

256  MB  each  node  ; 
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Figure  3:  Timings  for  a  ship  air  wake  case  on  several 
parallel  computers 


Figure  4:  Timings  for  a  ship  air  wake  case  on  SP2 


Results  and  Discussions 

In  the  helicopter /ship  interface  problem,  the  most  im¬ 
portant  data  for  coupling  the  airwake  solution  to  the  dy¬ 
namics  analysis  of  the  main  rotor  blades  of  helicopter  are 
the  mean  flow  field,  the  intensity  of  flow  perturbations, 
and  its  dominant  frequencies.  Such  an  approach  is  pre¬ 
sented  in  this  paper.  So  far  we  have  been  concentrating 
on  two  types  of  ships:  (1)  frigates  with  helicopter  land¬ 
ing  pads  on  the  deck  behind  the  hangar  and  (2)  aircraft 
carriers  and  LHA’s  with  several  helicopter  landing  spots 
on  the  deck  around  the  control  tower.  The  airwake  in¬ 
fluences  on  the  helicopter  are  quite  different  in  these  two 
cases.  For  frigates  the  flow  separation  area  behind  the 
hangar  cube  has  a  strong  effect  on  a  landing  helicopter, 
while  on  LHA’s  the  deck  leading  edge  vortex  and  sep¬ 
aration  are  the  key  flow  phenomena.  Fig.  5  is  some  of 
the  experimental  data  obtain  by  one  of  the  authors  (L. 
Long)  with  K.  Long  (PAX  River)  on  the  U.S.S  Saipan 
(an  LHA).  It  shows  the  flow  velocity  at  different  heights 
from  the  deck  in  the  center  plane.  This  was  obtained 
from  a  cup  anemometer  and  only  includes  the  effects  of 
velocity  in  the  longitudinal  and  normal  (to  deck)  direc¬ 
tions.  The  wind  was  at  36  knots  and  had  a  30  degree 
yaw  angle.  These  data  will  eventually  be  used  to  validate 
numerical  results  for  the  LHA  case. 


Figure  5:  Field  measurement  of  LHA 

In  this  paper  some  preliminary  simulations  have  been 
done  for  a  generic  ship  shape  (TTCP  ship).  In  fact,  this 
is  a  generic  frigate  model  and  it  is  shown  in  Fig.  6  in 
a  computational  mesh.  It  is  240  feet  long,  45  feet  wide 
and  55  feet  high.  It  was  chosen  because  there  will  be 
some  experimental  investigations  using  the  Scime  confi- 
gration.  It  is  acknowledged  that  the  ship  superstructure 
does  not  resemble  a  typical  frigate  superstructure  in  de¬ 
tail.  However,  from  an  aerodynamic  point  of  view  the 
airwake  should  still  be  representative  of  that  for  an  ac¬ 
tual  frigate  since  this  study  is  concerned  with  the  macro¬ 


scopic  flow  properties  and  large  scale  phenomena  in  the 
hangar  wake. 

Special  attention  is  given  to  the  helicopter  landing 
area  which  is  the  square,  aft  section  of  the  ship.  There  is 
a  20  feet  drop  down  to  the  landing  deck  from  the  hangar 
structure,  which  will  lead  to  vortex  shedding  over  the 
deck  causing  landing  approach  hazards. 

The  computational  grid  for  this  problem  is  201  x  109 
X  85  which  results  in  a  grid  resolution  of  two  feet  or 
less  in  each  direction  around  the  ship,  grid  stretching 
was  used  to  enlarge  the  domain.  So  far  both  mean  flow 
and  NLDE  simulations  were  based  on  the  same  grid  in 
order  to  avoid  three  dimensional  interpolation.  In  fact, 
the  NLDE  needs  much  fewer  grid  points  than  CFL3D. 

Mean  flow  simulation 

NASA  Langley  and  Ames  research  centers  have  de¬ 
voted  significant  resources  in  the  past  decades  to  devel¬ 
oping  modern  CFD  technology.  The  CFL3D  5.0  pack¬ 
age  from  NASA  Langley  is  used  here  to  simulate  the 
mean  flow  which  will  be  given  as  a  background  flow 
to  the  unsteady  flow  computation  of  NLDE.  The  code 
is  a  Reynolds- Averaged  thin-layer  Navier-Stokes  flow 
solver  for  structured  grids.  A  finite  volume  algorithm 
with  a  spatial-factored  diagonalised,  implicit  scheme  is 
used  in  discretization  of  the  partial  differential  equations. 
The  upwind-biased-differencing  using  the  flux-difference¬ 
splitting  technique  is  employed. 

From  the  experimental  results,  it  is  known  that  the 
flow  is  mostly  separated,  with  free  vortices  originating 
from  the  sharp  corners.  There  are  two  types  of  separa¬ 
tion:  one  due  to  viscosity  and  the  other  due  to  sharp 
corners  of  the  blocked  structures.  The  former  is  heavily 
influenced  by  the  Reynolds  number.  The  latter  is  purely 
an  inviscid  phenomenon,  independent  of  Reynolds  num¬ 
ber.  The  air  wake  is  greatly  influenced  by  both  of  them. 
In  this  mean  flow  simulation,  we  are  concerned  primarily 
with  the  inviscid  phenomenon  and  used  the  Euler  solver 
of  CFL3D.  The  TTCP  ship  computational  domain  is  di¬ 
vided  into  10  blocks. 

The  Mach  number  chosen  for  the  simulation  is  a  high 
wind  case.  The  incoming  flow  speed  is  41  knots.  The 
water  surface  is  assumed  to  be  a  hard  wall  boundary. 
Fig.  7  shows  the  contour  plot  of  velocity  magnitude  on 
the  surface  of  TTCP  ship.  This  is  a  zero  yaw  angle  case 
from  CFL3D  results.  The  asymmetry  property  of  zero 
yaw  angle  flow  is  captured  very  well.  The  flow  is  accel¬ 
erated  around  the  sharp  corners  and  there  are  several 
reverse  flow  regions  neeir  the  walls  close  to  each  corner. 
After  the  blocked  structures  there  is  massive  flow  sep¬ 
aration;  the  separation  line  is  clearly  shown  after  each 
block. 

Of  importance  to  the  landing  operation  is  the  flow 
condition  over  the  flight  deck.  Fig.  8  depicts  the  contour 


Figure  6:  Configuration  of  TTCP  ship  and  computational  mesh 


of  velocity  magnitude  at  the  ship’s  center  plane.  It  is 
shown  that  the  Icu-ge  region  of  recirculating  fiow  extends 
over  the  flight  deck  and  rises  higher  than  the  hangar. 
This  flow  region  is  in  the  landing  path. 


Figure  8:  Flow  speed  contour  on  center  plane  of  the 
TTCP  ship 

Fig.  9  shows  velocity  vectors  in  two  horizontal  planes 
4.75  feet  and  8.75  feet  above  the  flight  deck.  Our  numer¬ 
ical  results  are  compared  with  a  flow  pattern  obtained 
from  an  experimental  study^^  in  fig.  10.  It  shows  the 
flow  pattern  from  experiments,  where  four  distinct  flow 
regions  are  behind  the  hangar.  This  three  dimensional 
vortex  and  reverse  flow  has  very  low  speed  but  generally 


is  very  unsteady  and  yaw-dependent.  Comparing  to  the 
experiment,  the  physical  flow  features  are  well  captured 
by  the  simulation.  In  the  numerical  plots,  the  vortex 
pair  in  the  higher  plane  is  much  close  to  the  center  line 
and  the  hangar.  This  indicates  that  there  is  a  horse  shoe 
vortex  as  shown  by  the  topological  drawing  in  fig.  10. 

In  fig.  11  the  velocity  vector  on  the  deck  floor  is 
compared  with  fiow  visualization  results  for  the  TTCP 
ship.^^  There  are  differences  in  the  attachment  point 
and  the  position  of  vortex  center.  This  is  probably  due 
to  our  inviscid  appoach. 

Based  on  the  discussion  above,  in  fact,  the  flow  pat¬ 
tern  shown  represents  a  very  rough  mean  of  the  flow. 
They  are  intended  only  to  give  approximate  envelopes 
for  the  different  regions  and  provide  the  background  flow 
for  NLDE  simulations.  There  are  wild  fluctuations  about 
this  mean  flow.  The  flow  field  is  generally  very  unsteady. 
In  the  following  section,  the  results  from  NLDE  simula¬ 
tions  are  discussed  based  on  this  mean  fiow. 

Finally  in  fig.  12  the  performance  of  CFL3D  (a  serial 
code)  is  shown.  For  this  TTCP  ship  configuration  it 
takes  days  to  get  the  steady  state  results. 


Figure  7:  Flow  speed  contour  on  the  surface  of  TTCP  ship 


Figure  9:  Flow  velocity  vectors  at  two  horizontal  planes 
over  flight  deck 

Perturbation  simulation 

As  mentioned  before,  NLDE  has  been  developed  for 
solving  more  complex  geometries,  such  as  the  TTCP 
ship,  which  allow  multi-solid-boxes  inside  the  compu¬ 
tational  domain.  This  makes  the  boundary  condition 
implementation  difficult,  especially  combined  with  a  do¬ 
main  decomposition  parcdlel  technique.  To  overcome 
this  difficulty,  a  single  block  domain  was  chosen  for  the 
NLDE  code.  The  TTCP  ship  is  divided  into  88  solid 
boxes.  Characteristic  boundary  conditions  are  used  at 
the  surfaces,  edges  and  corners  of  these  boxes.  At  each 
time  step,  after  the  single  block  computation  is  finished, 
the  solid  box  wall  boundary  conditions  are  applied  to 
update  the  value  at  wall  grid  points. 


Figure  10:  Topological  drawings  from  experimented 
study 

The  high  wind  speeds  relevent  to  the  ship/helicopter 
interface  problem  arise  from  storm  centers  far  from  the 
actual  ship  and  are  called  neutrally  stratified.  This  wind 
condition  is  considered  at  our  inflow  boundary.  The 
principal  parameters  of  the  freest  ream  airflow  are  (1) 
the  mean  windspeed,  time-averaged  over  an  appropriate 
scale;  (2)  the  turbulence  intensity;  (3)  the  longitudinal 
(or  integral)  length  scale  of  the  turbulent  velocity  fluctu¬ 
ations.  Empirical  relationships  are  available  (ESDU  data 
items  74030,74031)  for  the  above  four  parameters  as  a 
function  of  the  mean  windspeed,  elevation  and  roughness 
length  scale. ^ 

Incoming  characteristics  with  source  terms  are  intro¬ 
duced  from  the  inflow  boundary  as  it  is  expressed  as  Li 
in  equation  (17)  -  (21).  The  magnitude  of  the  incoming 
disturbance  is  determined  by  the  turbulence  intensity. 
Its  spectrum  is  obtained  from  the  wind  spectrum  by  us¬ 
ing  a  random  walk(random  phase)  technique. 

Fig.  13  and  fig.  14  give  contour  plots  of  longitudinal 


Figure  11:  Predicted  surface  flow  velocities  compared  to  surface  oil  flow  images^^ 


Figure  12:  Timing  for  CFL3D  on  several  computers 


velocity  perturbations  u'  in  the  center  plane  at  two  dif¬ 
ferent  timesteps.  Fig.  15  gives  contour  plots  of  the  lon¬ 
gitudinal  velocity  intensity  From  this  unsteady 

result,  it  is  shown  that  large  perturbations  occur  around 
the  TTCP  ship  structure,  especially  in  the  area  after  the 
hangar  and  after  the  leading  edge.  In  the  field  far  from 
the  ship  the  flow  is  quite  steady.  Vortex  shedding  from 


the  hangar  can  be  clearly  observed. 

In  fig.  16  and  fig.  17,  contour  plots  of  vertical  and 
transverse  perturbation  velocity  are  shown.  In 

fig.  18  the  vertical  perturbation  intensity  \/w^  is  given 
in  the  center  plane.  High  instantanious  vertical  pertur¬ 
bations  are  found  in  the  region  just  after  the  hangar  trail¬ 
ing  edge.  Since  this  is  a  zero  yaw  angle  case,  the  trans¬ 
verse  perturbation  is  quite  weak  over  the  flight  deck.  By 
comparing  the  instantaneous  perturbations  in  fig.  13  and 
17,  the  perturbation  length  scale  of  longitudinal  pertur¬ 
bations  is  different  from  that  of  vertical  perturbations. 

The  unsteady  three-dimensional  flow  is  of  interest 
throughout  the  domain  but  in  particular  the  flow  un¬ 
steadiness  is  important  around  the  helicopter  landing 
deck.  Fig.  19  presents  a  contour  plot  of  perturbation 
intensity  in  a  horizontal  plane  17  feet  above  the  deck, 
where  the  helicopter  rotor  would  be. 

From  those  prelinimary  results,  the  unsteady  features 
of  TTCP  ship  air  wake  are  captured  qualitatively.  How¬ 
ever,  detailed  experimental  data  is  not  yet  available.  In 
the  meantime,  the  NLDE  code  is  being  improved  and 
prepared  for  quantitative  evaluation  and  analysis. 

Concluding  Remarks 

This  paper  presents  steady  and  unsteady  flow  field 


Figure  13:  Contour  of  instantaneous  longitudinal  velocity  perturbations  at  t 


Figure  14:  Contour  of  instantaneous  longitudinal  velocity  perturbations  at  i  =  t2 


predictions  for  frigate  class  ships.  A  nonlinear  distur¬ 
bance  equation  solver  has  been  developed  using  parallel 
computers.  The  parallel  performance  of  the  code  has 
been  compared  on  various  computers.  Our  present  re¬ 
sults  are  qualitatively  correct,  and  show  that  the  key  flow 
phenomena  can  be  captured  by  using  a  steady-state  code 
followed  by  the  NLDE  code.  Future  work  will  concen¬ 
trate  on  more  detailed  comparisons  to  experiment,  the 
inclusion  of  more  geometrical  features  of  the  ships,  and 
the  inclusion  viscous  effects. 
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